using System;
using System.Reflection;
using System.IO;
using NBody;
using NBody.Analysis;
using OptionParser;

[assembly: AssemblyTitle ("NBodyPotential")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Compute the potential.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class PotentialOptions : InputOutputGetOptions
{ 
    [Option ("-c", "Adjust for centre of mass")]
    public bool CentreOfMass;
     
    [Option ("-e", "Softening length")]
    public double Softening;
     
    [Option ("-f", "File for potential data")]
    public string PotFile;
    
    public PotentialOptions()
    {
        Softening = 0.0;
        PotFile = "";
    }
}

public class NBodyPotential
{
    public static void Main(string[] args)
    {
        PotentialOptions opts = new PotentialOptions();
        opts.ProcessArgs(args, ".rpot");
                        
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass(1e-3);
        
        if (opts.PotFile != "")
            Energy.SetPotentialFromData(ref s, DataFile.Table.Read(opts.PotFile));
        else
            Energy.CalculatePotentialShell(ref s, opts.Softening);
        
        using (StreamWriter sw = new StreamWriter(File.Open(opts.OutFile, FileMode.Create)))
        {
            int mod = (int)((double)s.NumParts / 10000.0);
            for (int i = 0; i < s.NumParts; i++)
            {
                if (i % mod == 0)
                {
                    EnergyParticle p = (EnergyParticle)s[i];
                    if (p.Potential < 0.0)
                        sw.WriteLine("{0} {1} {2}", Math.Log10(p.Radius), p.Potential, Math.Log10(-p.Potential));
                }
            }
        }
                
        /*int n = 100;
        DataFile.Data data = new DataFile.Data(n, 2);
        double rho0 = 0.00750241;
        double a = 62.9464;
        for (int i = 0; i < n; i++)
        {
            double r = Math.Pow(10.0, -1.0 + 4.0 * i / (n - 1.0));
            data[i, 0] = r;
            data[i, 1] = 4.0 * Math.PI * rho0 * a * a * a / r * Math.Log((a + r) / a);
        }
        
        data.Print("test.pot", new DataFile.Transform(Math.Log10));*/
    }
}
